**objetivo de aprendizaje” Este caso práctico muestra como leer y graficar datos espaciales en R. BLA BLA

Tarea 1: Abrimos Rstudio

Tarea 2: Leemos (instalamos si no los tenemos) los paquetes requeridos

# rm(list = ls()) #limpiamos la memoría

# library(climaemet) # meteorological data
library(mapSpain) # base maps of Spain
library(classInt) # classification
library(raster) # raster handling
library(sf) # spatial shape handling
library(sp) # spatial shape handling
library(gstat) # for spatial interpolation
library(geoR) # for spatial analysis
library(dplyr) # data handling
library(ggplot2) # for plots
# library(tidyverse) # collection of R packages designed for data science

Tarea 3: Describimos los datos objeto de estudio

El conjunto de datos contiene el nivel de temperatura del aire en España el 8 al 13 de enero de 2021.

Los datos han sido descargados del paquete climaemet en R. climaemet permite descargar los datos climáticos de la Agencia Española de Meteorología (AEMET) directamente utilizando su API (para más detalle ver citar articulo).

This dataset has several dataframes:

Tarea 4: Lectura de los datos

mydata <- read.csv("my_data_tmin_df.csv")
mydata <- mydata %>% select(X:date) # quitamos la primera columna que tiene id

Q1: ¿Qué información tiene “mydata?”

head(mydata)[1:3, ] # muestra las tres primeras filas
#>           X        Y tmin       date
#> 1 -2.956389 35.27639 13.6 2021-01-08
#> 2 -5.346944 35.88861 11.6 2021-01-08
#> 3 -5.598889 36.01389 11.1 2021-01-08
tail(mydata) # muestra las seis últimas filas. Ver fecha de los días.
#>              X        Y  tmin       date
#> 1291 -3.528333 39.49194 -13.4 2021-01-13
#> 1292 -3.742778 41.66583 -13.9 2021-01-13
#> 1293 -1.410000 41.11444 -14.9 2021-01-13
#> 1294 -1.124167 40.35056 -15.5 2021-01-13
#> 1295 -1.293333 40.92611 -18.0 2021-01-13
#> 1296 -1.878889 40.84167 -19.9 2021-01-13

Tarea X. Pasar mydata a un objeto de la clase espacial geoR

mydata.geoR <- mydata %>%
  filter(date == "2021-01-08") %>%
  as.geodata(
    coords.col = 1:2,
    data.col = 3
  )
#> as.geodata: 4 replicated data locations found. 
#>  Consider using jitterDupCoords() for jittering replicated locations. 
#> WARNING: there are data at coincident or very closed locations, some of the geoR's functions may not work.
#>  Use function dup.coords() to locate duplicated coordinates.
#>  Consider using jitterDupCoords() for jittering replicated locations
summary(mydata.geoR)
#> Number of data points: 215 
#> 
#> Coordinates summary
#>             X        Y
#> min -9.291389 35.27639
#> max  4.215556 43.78611
#> 
#> Distance summary
#>      min      max 
#>  0.00000 13.85144 
#> 
#> Data summary
#>        Min.     1st Qu.      Median        Mean     3rd Qu.        Max. 
#> -14.9000000  -4.5500000  -0.5000000  -0.6060465   3.5000000  13.6000000 
#> 
#> Duplicated Coordinates
#>     dup          X        Y data
#> 40    1 -0.3663889 39.48056  5.0
#> 46    2 -0.3663889 39.48056  4.3
#> 139   3 -1.7922222 43.35694 -2.9
#> 142   4 -1.7922222 43.35694 -3.0
#> 67    5 -3.8005556 43.49111  2.5
#> 68    6 -3.8005556 43.49111  2.4
#> 130   7 -5.8741667 43.35333 -2.1
#> 131   8 -5.8741667 43.35333 -2.1
plot(mydata.geoR)

DIEGO, y yo por qué se que España es esto 4326…

Tarea X. Pasar mydata a un objeto de la clase espacial sf

mydata.gstat <- st_as_sf(mydata, coords = c("X", "Y"), crs = 4326)
summary(mydata.gstat)
#>       tmin             date                    geometry   
#>  Min.   :-25.200   Length:1296        POINT        :1296  
#>  1st Qu.: -4.500   Class :character   epsg:4326    :   0  
#>  Median : -0.800   Mode  :character   +proj=long...:   0  
#>  Mean   : -1.144                                          
#>  3rd Qu.:  2.600                                          
#>  Max.   : 13.600

Tarea X. Decido un formato para el análisis. gstat

mydata <- mydata.gstat

Q2: ¿Dónde tengo medidos los niveles de la temperatura mínima?

plot(mydata$geometry, pch = "+")

Tarea 5: Dibujemos las estaciones de monitoreo de la temperarua mínima en un mapa de España. Ámbito espacial.

Obtenemos el contorno de España (quitamos las Islas Canarias por simplicidad) blaalblblbl

library(mapSpain)
# sf object
ESP <- esp_get_ccaa(epsg = 4326) %>%
  filter(ine.ccaa.name != "Canarias") %>%
  st_union()
plot(ESP) # Dibujamo el mapa de España menos las Islas Canarias

Q3: ¿Tengo el Sistema de referencia de coordenadas (CRS) de las estaciones de monitoreo en la misma proyección que el contorno de España?

st_crs(mydata) == st_crs(ESP)
#> [1] TRUE

Dibujamos las estaciones de monitoreo con el contorno de España

ggplot(ESP) +
  geom_sf() +
  geom_sf(data = mydata) +
  theme_light() +
  labs(
    title = "Estaciones de monitorieso AEMET en  España",
    subtitle = "excluyendo las Islas Canarias"
  ) +
  theme(
    plot.title = element_text(
      size = 12,
      face = "bold"
    ),
    plot.subtitle = element_text(
      size = 8,
      face = "italic"
    )
  )

Q4. Mis datos y el contorno de España están en el mismo CRS, pero ¿es el adecuado?

DIEGO, POR QUÉ TRANFSORMAMOS, PARA PASAR A UTM E INTERPRETAR EN METROS?

mydata <- st_transform(mydata, 25830)
ESP_utm <- st_transform(ESP, 25830)

Tarea 6: Representamos la variable temperatura mínima tmin para el día 8 de enero de 2021.

mydata_8enero <- mydata %>%
  filter(date == "2021-01-08") # seleccionamos el día y la variable
dim(mydata)
#> [1] 1296    3

# plot(mydata_8enero) #opcional
plot(mydata_8enero["tmin"], main = "Temperatura mínima (8-enero-2021)", pch = 8)

Podemos utilizar el ámbito espacial, el contorno de España para graficar y contar la historia de la Filomena un poco mejor.

Q5. El mapa ha quedado muy claro. Vemos como los datos nos cuentan la historia de Filomena en aquellos sitios donde se tomaron mediciones, pero ¿podríamos tener un mapa de interpolación para tener una estimación de la temperatura mínima en las partes donde la AEMET no tiene estación de monitoreo?

Tal y como se avanzó en teoría, parece lógico pensar que aquellos puntos que estén cerca tendrán valores similares así que tomemos ventaja de la dependencia espacial y utilicemos un método determinista, como la Distancia Inversa Ponderada, comúnmente conocido por su acrónimo inglés IDW (Inverse distance weighted), el cual es uno de los métodos más simples para llevar para llevar a cabo una interpolación espacial.

En primer lugar necesitamos definir la superficie (en forma de malla) donde queremos interpolar

# Creamos una malla de 5*5 km (25 km2)
grd_sf <- st_bbox(ESP_utm) %>%
  st_as_sfc() %>%
  st_make_grid(
    cellsize = 5000,
    what = "centers"
  )

# Convertimos a un objeto sp object
grd <- as(grd_sf, "Spatial") %>%
  as("SpatialPixels")

Graficamos la superficie para ver exactamente qué hemos construido

ggplot(ESP_utm) +
  geom_sf() +
  geom_sf(data = grd_sf, size = 0.01, col = "red", alpha = 0.5) +
  geom_sf(
    data = mydata_8enero, # select(indicativo) %>% unique(),
    aes(fill = "AEMET Stations"), size = 5, shape = 21,
    color = "blue"
  ) +
  scale_fill_manual(values = adjustcolor("blue", alpha.f = 0.2)) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  labs(
    title = "Cuadrícula espacial para interpolar",
    fill = ""
  )

Cambiamos el conjunto de datos a la clase sp para el análisis espacial

mydata_8enero_sp <- as(mydata_8enero, "Spatial")
class(mydata_8enero_sp)
#> [1] "SpatialPointsDataFrame"
#> attr(,"package")
#> [1] "sp"

LLevamos a cabo la interpolación

tmin_idw <- gstat::idw(tmin ~ 1, # Indicamos la variable que queremos interpolar
  mydata_8enero_sp, # Indicamos el conjunto de datos donde está la variable
  newdata = grd, # Indicamos los puntos a interpolar
  idp = 2.0 # Especifica la potencia de la IDW
)
#> [inverse distance weighted interpolation]

DIEGO, mirame MASK, no entiendo el objeto que hay que poner

Representamos la interpolación con un mapa y mapa de contorno muy utliazado para representar datos espaciales

# To raster and mask the grid to the shape of Spain
tmin_idw_raster <- raster(tmin_idw) # %>% mask( ) ##ERROR
plot(tmin_idw_raster)
contour(tmin_idw_raster, add = TRUE)

Representamos en 3D el mapa anterior, muy utilizado también en datos espaciales

DIEGO, se útiliza mucho en espacial pero no se si este caso es muy ilustrativo….r

persp(tmin_idw_raster)

0.1 Para las extensiones

Este vale para st